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We apply the adaptive time-dependent Density Matrix Renormalization Group method (tDMRG) 
to the study of transport properties of quantum-dot systems connected to metallic leads. Finite-size 
effects make the usual tDMRG description of the Kondo regime a numerically demanding task. We 
show that such effects can be attenuated by describing the leads by "Wilson chains" , in which the 
hopping matrix elements decay exponentially away from the impurity {t A""/"). For a given 
system size and in the linear response regime, results for A > 1 show several improvements over the 
undamped, A = 1 case: perfect conductance is obtained deeper in the strongly interacting regime 
and current plateaus remain well defined for longer time scales. Similar improvements were obtained 
in the finite-bias regime up to bias voltages of the order of the Kondo temperature. These results 
show that, with the proposed modification, the tDMRG characterization of Kondo correlations in 
the transport properties can be substantially improved, while it turns out to be sufficient to work 
with much smaller system sizes. We discuss the numerical cost of this approach with respect to the 
necessary system sizes and the entanglement growth during the time-evolution. 

PACS numbers: 73.63.Kv,73.23.-b, 72.10.Fk, 72.15.Qm 



I. INTRODUCTION 



The current excitement in the condensed matter and 
materials science community surrounding the study of 
nanoscale transport stems from both the potential appli- 
cability in molecular electronic devices^ and the possibil- 
ity of designing nanostructures to realize quantum impu- 
rity Hamiltonians. Hallmark experimental achievements 
include the observation of the Kondo effect in quantum 
dotsr'i^ molecules^ nanotubeS)^ and non-Fermi liquid 
behavior in quantum dot structuresii 

While transport is an intrinsic non-equilibrium sit- 
uation, in the linear response regime transport coeffi- 
cients are commonly derived from equilibrium correla- 
tion functions i^i^ A prominent numerical tool for describ- 
ing the equilibrium Kondo regime is Wilson's numeri- 
cal renormalization-group method (NRG). In the orig- 
inal NRG formulation for the Kondo modelji^iii Wil- 
son showed that the contribution from band states ex- 
ponentially close to the Fermi energy needs to be taken 
into account in order to capture the correct properties 
of the ground state. For this reason, standard tight- 
binding numerical approaches face a formidable challenge 
in addressing this problem: finite-size effects set a min- 
imum energy scale, the level spacing, below which the 
calculation cannot capture the crossover to the Kondo 



Wilson proposed a combination of two elements to 
handle this problem: (i) A discretization procedure of 
the metallic band, leading to a mapping into a impu- 
rity connected to a one-dimensional tight-binding chain 
with exponentially decaying hoppings. We will refer 
to such leads as Wilson chains in this work, (ii) A 
non-perturbative renormalization procedure that probes 
successive energy scales by recursively diagonalizing the 
Hamiltonian and keeping the relevant states at each scale. 

Recent theoreticali^^'i^^i20^'22^^^^^^^^^ 
and experimental efforts^ aim at observing and mod- 
eling genuine non-equilibrium physics. A particularly 
important question is under what conditions steady- 
state situations can be reached in numerical simulations, 
and promising results have been obtained using time- 
dependent approachesi ^^i^^'^^i^^'^'^i^^ Such ideas have 
been pursued using both the density matrix renor- 
malization group (DMRG) technique^i^ii^i^^ and the 
NRG,— i^i^^i^ in the former case utilizing the adaptive 
time-dependent DMRG^i^^ 

Moreover, the incorporation of ingredients of DMRG 
into NRG and vice versa has led to a significant extension 
of both methodsi ^^i^^i'^^i^-'^i'^^i^'^ As a prominent example, 
this includes the use of Wilson chains in DMRG for the 
description of the Kondo regime of the Anderson impu- 
rity model in Refs. [40ll4ll . where in Ref. [4l| the common 



2 



mathematical structure of NRG and (single-site) DMRG 
in terms of matrix-product states has been exploited. 
Recently, a similar idea has been successfully explored 
within a cluster-embedding approach, resulting in the 
development of the so-called logarithmic discretization 
embedded cluster approximation (LDECA).— 

The advantage of DMRG is its flexibility: it is in princi- 
ple possible to model complex interacting regions^ or to 
incorporate interactions into the leads (see, e.g., Ref . [ish . 
Moreover, it is the numerical method of choice for one- 
dimensional bulk systems, and it allows for the calcula- 
tion of extended correlation functions in a straightfor- 
ward way. For the description of transport phenomena, 
there is no restriction to work in the small bias regime, 
as finite biases can be incorporated into time-dependent 
simulations 

In transport investigations based on DMRG, sev- 
eral groups have introduced modifications in the con- 
tact leads, such as the logarithmic discretization, to 
improve the results of either the ground state^>^ or 
tDMRG calculations. These also include damped bound- 
ary conditions22ii& and a momentum-space representa- 
tion of the leads4i In the latter work, an interacting res- 
onant level model has been studied with tDMRG and, 
effectively, a logarithmic discretization of the leads has 
been used. Working with different kind of leads, while 
preserving the main physical properties of the system, 
has thus proven to be a promising path that we will fur- 
ther pursue in this work. 

It is the purpose of this paper to perform a numerical 
real-time analysis of Kondo correlations in quantum dot 
problems using tDMRG and Wilson chains. We show 
that, compared to a previous study by some of us^^ a 
correct description of transport through a quantum dot 
can be obtained deeper into the Kondo regime, and using 
much smaller system sizes. While a tDMRG analysis of 
the Kondo regime based on a real-space description is 
hampered by finite-size effects in the leads, we show that 
an appropriate choice of hopping amplitudes in the leads 
nicely circumvents such problems. 

In this sense, the discretization scheme proposed by 
Wilson is a natural choice: the noninteracting tight- 
binding chain becomes an effectively metallic system 
with reduced level spacings, corresponding to a subset 
of states directly coupled to the impurity in the contin- 
uum limit i^ii^ In addition, it turns out to offer advan- 
tages in the time-dependent description as well: for a 
given system size, it substantially increases the charac- 
teristic time scales over which a constant current flows 
between the leads, before reaching the system's bound- 
ary. This turns out to be crucial deeper into the Kondo 
regime (i.e., for small Kondo temperatures), as constant 
currents sustained over longer time scales are necessary 
to access the regime of coherent transport through the 
sharp, resonant Kondo state.— i^i^ We also discuss and 
study the influence of the discretization at finite-bias. In 
the low-bias regime, the logarithmic discretization gives 
the correct description of the Kondo regime, at large bias. 



(where the Kondo effect is effectively destroyed), a lin- 
ear discretization of the leads and larger chains should 
be used. 

The paper is organized as follows: In Sec. [IT] we in- 
troduce our approach, using the example of a single 
quantum-dot connected to metallic leads out of equi- 
librium, treated with the tDMRG method. Results for 
the time-dependent currents and charge transfer are dis- 
cussed in Sees. |lll] and IIVI respectively. Particular em- 
phasis is devoted to describing the improvements ob- 
tained by choosing a A > 1 model over the A = 1 case, 
and to the dependence of our results on both the dis- 
cretization parameter and the system size. Moreover, 
we discuss the conductance results, taken from the time- 
dependent data, in Sec. |Vl In Sec. IVIl we provide a 
discussion of how to approach the finite-bias regime by 
combining tDMRG data obtained from discretized and 
tight-binding leads. We present a summary in Sec. IVIII 
Two appendices conclude this work: App.EI contains our 
results for the noninteracting case, and App.lBl provides a 
discussion on the computational aspects of our approach 
and the entanglement growth during the time-evolution. 



II. MODEL AND SET-UP 

As a case study of the proposed modification of the 
method, we apply tDMRG to a model representing a 
single quantum dot connected to metallic leads. The 
equilibrium and linear-response properties of this system 
are well known and provide a natural benchmark against 
which we can compare the tDMRG results. 

The quantum dot is modeled by a Hubbard site with 
an on-site interaction U and a gate potential Vg coupled 
to noninteracting tight-binding chains, representing the 
leads, as depicted in Fig. [1] The dot-lead couplings are 
labeled by t' . The full Hamiltonian reads H = i/Dot + 

-f^Dot-Lcads + -f^Lcads , with 
^^Dot-Loads = -t' ^ {cl^Ci^a,a + h.C.^ , 

<J,a—R,L n—1 

(1) 

where cjj q, creates an electron with spin a at site n in 
lead a {n = is the dot site), Na counts the number 
of sites in lead a, and •fiia = c\^c^^; hi = fii-^ + hi^ 
with i = 0, {«,«}. The total number of sites is thus 
= Nl + 1 + Nr. Our results were obtained using 
"even-l-odd" chains, with Nl = N/2, Nr = N/2 - 1. 
Note that finite-size effects due to different cluster types, 
as discussed in Ref. H^, vanish at sufficiently large A. 

In the spirit of Wilson's discretization scheme, we con- 
sider the hopping matrix elements t„ in the leads to decay 
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FIG. 1: (color online) Schematic representation of the model 
describing a quantum dot (large circle, blue in the online ver- 
sion) connected to left and right leads. 

exponentially as 

t„=toA-"/2, (2) 

where A > 1 is the "discretization parameter" from Wil- 
son's original formulation. 

Unless otherwise noted, in the following we set to = 1, 
U = 1. Vg ~ ~U/2 (particlc-holc symmetric point). We 
focus on results for = 32 sites for A > 1 and up to 
N = 128 sites for A = 1. We usually work at half-filling of 
the whole system. In the range of parameters considered, 
the equilibrium ground-state properties for A > 1 do not 
change significantly upon increasing N beyond a certain 
N*{K) < 50 at a given A, as we have numerically verified 
for a few cases. This is consistent with NRG runs for an 
Anderson model with similar parameters, which reach 
the Kondo fixed point with less than 50 iterations for 
A -2- 3. 

Details on the tDMRG can be found in Refs. [l^IH. 
We use a Trotter-Suzuki breakup of the time-evolution 
operator and typical time steps of 5t = 0.01-0.1. A 
larger time-step St ~ 0.4 is sufficient when deeper in the 
Kondo regime (large U/t'^) as resonant transport is dom- 
inated by a small energy scale, the Kondo temperature, 
corresponding to long time scalesJ^i^ Note that we de- 
note time by the symbol r and it is measured in units of 
h/tQ. The truncated weight 6p during the time-evolution 
is typically kept below 10"'' (see App.|B]for a discussion). 

In order to drive a current, we first compute the ground 
state of the system without a bias. The bias is applied 
as 

-"bias = 2^ nji,n — nL,n (3) 

n— 1 n— 1 

in the leads at time r and we then evolve in time under 
the dynamics oi H + i/bias- We typically work at a small 
bias of AV ~ 0.005 (the finite-bias case is discussed in 
Sec. IVI[) . The current J(t) is measured as the average 
over the expectation values of the local current operator 

= jt'^(4_o.ci,Q,<T - h.c), (4) 

cr 

on the links connecting the dot to the leads. This means 
that we first take its expectation value in the time- 
dependent wave-function and then average over the two 
local currents on the links directly connected to the dot. 
We have tried other spatial forms for the applied bias, 
e.g., a broadened step functioni^ This mostly affects the 
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FIG. 2: (color online) Conductance J(r)/AV for U = 1, t' = 
0.4, Vg = -17/2. (a) Fixed A(= 1), and = 16, 32, 64 and 
128; (b),(c): Fixed A^. (b) A = 16 and (c) A = 32 with 
A = 1, 2^/*, V2, and 2. 



short-time transient behavior, but leaves unaffected the 
average value of the current taken over time (see Scc.|Vl). 
Further details on the set-up can be found in Ref. [20. 



III. TIME-DEPENDENT CURRENTS 

Figure [2] shows the current J(r) (in units of e'^/h) 
as a function of time and divided by the external bias 
AV, at f/ = 1, t' = 0.4. This corresponds to a ratio of 
U/r = 3.125, where F = 7rpieads(^')^ is the hybridization 
parameter and pieads = 2/(7rto) is the density of states of 
the leads in the limit of A = 1 and long chains. 

Figure[2i;a) contains the results for N = 16, 32, 64, 128 
at A = 1, reproducing those of Ref. The other two 
panels display J{t)/AV for (b) = 16 and (c) A^ = 32, 
computed with A = 1, 2'^^^, ^2, and 2. 

The comparison between Fig.l^Ja) and (c) is revealing. 
In the A = 1 case, the conductance plateaus become 
longer and higher as the system size increases^S with the 
"plateau length" (= Tgt) increasing linearly with N. For 
A^ = 128 sites, a nearly perfect conductance plateau with 
G(r) w Go is obtained for these parameters. A finite-size 
scaling analysis, done in Ref. [20 for U/T = 3.125, shows 
that G Go for A^-^ 0. 

As previously argued in Ref. [13, this plateau signals 
the formation of a Kondo state in the system formed by 
the quantum dot and the leads. A similar Kondo conduc- 
tance plateau can be obtained by increasing A while keep- 
ing the system size constant, as shown in Figs. [Ifb) and 



4 



(c). Interestingly, for N = 32, a plateau with an average 
conductance of G « Gq can be obtained with A = 2. 
Thus, an accurate description of the Kondo regime can 
be obtained using a relatively small system by taking 
A> 1. 

Notice that the J{t)/AV curve for A = 2 and iV = 32 
is similar to the corresponding one for iV = 16: an in- 
crease in system size had little effect in the current in this 
case, as opposed to the A = 1 curves. This weak depen- 
dence of J(t) on TV for large enough N and A is a general 
feature of the A > 1 case, and it is a consequence of the 
exponential decrease of the hopping matrix elements tn 
at large n. 

This allows us to formulate some key results that can 
be inferred from this plot: (i) By comparison of Fig. [2] (a) 
and Fig.[2](c), we find that with discretized leads, a four 
times smaller system size is sufficient to obtain roughly 
the same average current, (ii) At A > 1 and for a given 
the current plateaus are longer in time and, on the 
time-scales simulated, we do not observe a recurrence 
(bouncing current) in the case of TV = 32. (iii) Oscil- 
lations about the average value of the current tend to 
increase with A. 

The small time-scale oscillations are a general feature 
of the A > 1 case, which emerges in the U = case as 
well, as we have verified with exact diagonalization (see 
App. We have conducted extensive checks to rule 
out either the truncation error during the time evolution 
or the size of the time-step as possible sources of such 
oscillations at finite U. In fact, the position of the peaks 
and valleys of the oscillations as seen in Figs. [DJb) and 
(c) are A-dependcnt: for a given A, the position of these 
peaks and valleys remains practically constant even when 
other parameters, such as t' or Vg, are modified, while it 
changes as A is varied (see App. [A| . 

The oscillations are reminiscent of the so called "cur- 
rent ringing" in mesoscopic transport Similar current 
ringing effects have been observed in previous tDMRG 
studies of noninteracting systems away from half-filling 
or at finite biasi^'2^'^ In the present case of both a finite 
U and A > 1, the oscillations are present even at half- 
filling and become more prominent at larger values of A. 
We stress that the oscillations seen in the A > 1 case are 
an artefact of the discretization and are thus not a prop- 
erty of the A 1, N ^ oo limit. In practice, the oscilla- 
tions can be reduced by using the so-called z-trick,iiii^i^ 
which we illustrate in App. El for the noninteracting case. 

We attribute the increase of the average current on 
small chains at larger A to a combination of mainly two 
factors: (i) an effective reduction in the mean-level spac- 
ing in the leads in equilibrium and, quite importantly, (ii) 
an increase in the duration of the non-equilibrium cur- 
rent plateaus above a characteristic Kondo time scale. 
As the exponential decay of hoppings leads to an expo- 
nential decrease in the velocity at which charges move 
far away from the dot, those charges get trapped at the 
leads and a recurrence, i.e., a reversal of the current's 
sign, is not observed. We will elaborate on this point in 



0.05 




time 

FIG. 3: (color online) Charge transfer for A = 1, iV = 32 
(dashed lines), A = 1, A'' = 128 (open circles) and A — \/2, 
A = 32 (solid line). The parameters are U = 1, Vg = —U/2, 
and t' = 0.4 {U/T — 3.125). The maximum in the charge 
transfer for N = 32, A = 1 indicates a reversal in the current. 

the following section. 

Before turning to the calculation of conductances from 
our time-dependent data, we will discuss the charge pro- 
files and charge transfer during the time-evolution. 

IV. CHARGE PROFILE AND CHARGE 
TRANSFER 

Since no dissipative terms are included in the Hamil- 
tonian ([1]), the total charge is conserved at all times. 
Thus, the existence of a net current signals the transfer 
of charge from one lead to the other. As time progresses, 
a saturation point might be reached, opening the possi- 
bility for the current to decrease and reverse sign, and to 
transfer the excess charge back to the original lead.— i^i 

This mechanism is shown, for instance, in Fig. [2Ka) 
(A = 1): for iV = 16 the current reverses sign aroimd 
r = 9 while for = 32 the sign reversal occurs at r = 18. 
Such sign reversal also occurs for N = 128 and A = 1, at 
times T ^ 25. 

Figure [3] shows the charge transfer into the left lead 
(l), defined as (^^(r) = (n^r))): 

AnL(T) =^[n,(T)-n,;(G)] . (5) 

Notice that this quantity is related to the time- 
integrated current through the dot: it hence reaches a 
maximum whenever the current changes sign. For A = 1 
and N ~ 32, An;(T) reaches a maximum (ArT-f^max ~ 
0.01) at r = 18 while for A = \/2, it increases to about 
four times that value. Remarkably, in order to obtain a 
similar charge transfer enhancement with "regular" leads 
(A = 1), a four-fold increase in the system size is needed 
(open circles in Fig. [3]). 

The approximately linear increase in An;(r) is corre- 
lated with the plateau in the current in both cases (com- 
pare with Fig. [2]). We identify this as a general feature 
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FIG. 4: (color online) Charge profile of the 16-1-15 chain 
(single dot), Vg = -U /2: (a) A = 1 and (b) A = ^/2. The dot 
site is indicated by arrows in both plots. 
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FIG. 5: (color online) Scaling with 1/A of the conductance 
G for (a) N =ie and (b) = 32, Vg = -(7/2, and different 
values of fZ/F. 



panel): a larger charge is transferred to the left lead and 
no reflux is noticed. More importantly, the charge tends 
to accumulate toward the edge of the leads, with strong 
Friedel oscillations. This can be intuitively understood 
from the exponential decrease in the couplings: the weak 
coupling of the end sites with the remaining of the chain 
makes them good "charge traps" . 



of the introduction of the decaying hopping matrix el- 
ements for a given system size N: an enhanced charge 
transfer over longer time scales. 

This indicates that the exponential decay in the hop- 
pings increases the maximum charge that can be "stored" 
in the leads, or, in other words, it provides an increase 
in their effective "capacitance" . As a consequence, even 
small systems can hold a larger amount of charge without 
reversing the current, leading to longer constant-current 
plateaus. 

The effect of A > 1 can be illustrated by considering 
the time evolution of the charge profile. Figure |4] shows 
the charge n; (t) on each site I of the chain plotted against 
time for a chain of = 32 sites. The top panel shows 
the A = 1 case: charge is initially transferred from the 
left to the right lead and back, leading to an oscillation in 
ni(r) with a maximum around r = 18, as expected. The 
charge versus time diagram clearly shows a light-cone 
with a wave-front that propagates at the Fermi velocity. 
Notice that the excess charge is mostly accumulated in 
the vicinity of the dot, which we expect to modify the 
leads' density of states "seen" by the dot. 

This is in sharp contrast with the A = V2 case (bottom 



V. THE CONDUCTANCE 

We now turn to the linear conductance G(t) = 
J(t)/AF, expressed in units of Go = 2e^//i. We will dis- 
cuss the dependence of G(t) on A and N, at the particle- 
hole symmetric point. For this purpose, we refer the 
reader back to Fig. [2l Taking the example of A = 2, we 
see that no significant difference exists between the con- 
ductance curves for A^ = 16 and A^ = 32. In general, for 
a given A there exists a certain system size above which 
no significant changes in properties, such as the ground 
state energy or the conductance, take place as extra sites 
are added to the system. 

In order to study the dependence of the conductance 
with other parameters, we calculate the average conduc- 
tance G = {G(t))t over a plateau, e.g., those displayed 
in Fig. [21 This procedure carries an intrinsic uncertainty 
which depends on the truncated weight in the DMRG 
time evolution and, more importantly, on the dispersion 
of G(r) around the average due to the current ringing 
effects at larger A. While the former can be constrained 
below a target value by increasing the number of states 
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FIG. 6: (color online) Average conductance G = (G)t vs. 
U /T for A = 1, 2 and 3, and N=32 (dashed line is a guide to 
the eye). For A > 1, the averages are taken over time intervals 
"^li r < Tmax. For U/T <7, maximum simulation times are 
Tmax ~ tk + 100. For (//F > 8 (included for illustration), 
tk > Tmax, heucc G < Gq. luscti Kondo time tk (calculated 
with NRG) vs. U/V. See the text for details. 

that are kept during the time evolution (see the discus- 
sion in App. |B|, the latter is intrinsic for A > 1. We 
estimate such uncertainty by computing 

5G = v/(G(r)2),-(G(r))2, (6) 

and indicate it as error bars in the figures. We remark 
that the main contribution to 5G in the plots comes 
from the current-ringing oscillations (see the analysis in 
App.EJ. 

Figures EJa) and (b) show the scaling of G = (G(T))r 
with 1/A for different values of U/T, and for = 16 
and 32, respectively. The scaling is more conclusive for 
TV = 32: G ^ Go as 1/A decreases, for U/T < 5. Most 
importantly, Fig. EJb) establishes the convergence of the 
conductance in 1/A (obtained at a fixed system size) to 
the correct result, namely perfect conductance. 

Figure [S] depicts results for G/Gq as a function of U /F, 
for N = 32, = 0.005, and A = 1, 2, and 3. With 
A > 1, wc obtain perfect conductance up to [//F sa 7. 
This constitutes a considerable improvement over the 
A = 1 case with = 32 (also shown in Fig. (5]), for 
which perfect conductance plateaus were not observed for 
nonzero C//F. Furthermore, we stress that the A > 1 ap- 
proach also gives more well defined plateaus of constant 
currents, which in practice makes easier the averaging of 
J(t)/AF over time. 

The improvement previously discussed is anchored on 



a combination of two key elements in the A > 1 case: (i) 
an effective reduction in the level spacing of the metal- 
lic leads in equilibrium^ and (ii) the suppression of the 
current reversal, which is a consequence of the reduced 
velocity in the leads when A > 1. Both points are related 
to the exponential decrease in the chain hoppings, which, 
in Wilson's scheme, can be traced back to a representa- 
tion of the continuum of states directly connected to the 
quantum dot. 

Point (ii) is important for the following reason: in res- 
onant transport, the typical time scale for reaching the 
steady state is inversely proportional to the width of the 
resonance The typical time scale associated with the 
development of the perfect conductance plateaus associ- 
ated with the Kondo state is thus tk = h/Tx^^-^ It is 
then crucial that the current does not reverse sign before 
times of order t/<- have been reached. 

As explained in Sec. IIIIl for A = 1 the plateaus last 
over time intervals Tgt cx; N and a compromise must be 
obtained between Tk and N such that the condition Tgt > 
Tk is fulfilled. For A > 1, this condition can be met by 
increasing A (instead of N'), since Tgt increases with A, 
as discussed in Sec. IIIIl Constant currents corresponding 
to perfect conductance can, in principle, be reached for 
A values large enough so that Tst(A) > tk- In this sense, 
this requirement marks the regime for which the steady 
state of this problem can be numerically simulated. 

Once the condition Tst(A) > tk is fulfilled, one still 
needs to run the tDMRG algorithm over time scales of 
order tk to obtain the Kondo plateau. Thus, for higher 
values of U/T (i.e., higher tk)-, calculations over longer 
time scales are necessary in order to reach nearly perfect 
conductance plateaus in J(t).— This is, due to the entan- 
glement growth in a global quench^^i^ the true limita- 
tion of the method for Kondo problems. Fortunately, the 
entanglement growth turns out to be softer at large val- 
ues of U/T (see App. [B]), enabhng us to observe G ~ Gq 
for up to U/T ~ 6 and relatively short [N = 32) chains. 

This argument is further supported by quantitative es- 
timates for Tk obtained with NRG (inset in Fig.lH]). We 
performed NRG calculations for the Anderson model and 
extracted the Kondo temperature (and thus tk) from 
the magnetic susceptibility curves for different values of 
U /Fj^ For the parameters in Fig. [SI we obtain tk < Tgt 
in the regime where nearly perfect conductance is seen 
in the tDMRG curves {tk ^ 16 for U/T = 3.125 and 
Tk — 55 for U/T = 5, in units of h/to)^ For higher val- 
ues of U /F, Tk becomes exponentially large. In particu- 
lar, for C//F ~ 7, Tk calculated from NRG becomes of the 
order of the maximum time scales used in our tDMRG 
simulations. This explains the noticeable deviation of G 
from the Kondo value for U/T > 7. 

In short, the tDMRG results obtained with A > 1 con- 
stitute a considerable improvement over the A = 1 case, 
as shown in Fig. [6] for N = 32. This plot illustrates the 
range of parameters for which the Kondo regime is ac- 
cessible with tDMRG and A > 1, as well as the typical 
system sizes. 
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VI. FINITE BIAS 

In this Section, we address the case of a current 
through a quantum dot in the Kondo regime driven by 
a finite bias. Although a comprehensive theoretical un- 
derstanding of this nonequilibrium regime is yet to be 
achieved, one commonly expects the applied bias AV 
to disrupt the Kondo state for bias voltages larger than 
Tx while Kondo-like properties are only marginally af- 
fected for AV < Tif^^i^ai^aiSi We investigate the trans- 
port properties of the system in these two regimes. In 
the following; we fix the parameters to U/T = 3.125 
at the particlc-holc symmetric point {Vg = —U/2), for 
which wc have independently determined Tr- from NRG 
calculations^ (see the inset in Fig. [6]). 

Qualitatively, we expect the linear regime (i.e., nearly 
perfect conductance) to extend up to biases AV ^ Tk- 
The current versus bias curve should then smoothly drop 
to zero at larger biases. We argue that, in the spirit of the 
previous sections, in the regime AV < Tk, the finite-bias 
regime should best be explored with discretized leads, for 
which the Kondo state is better described. 

By contrast, in the opposite limit of AV 3> Tk, the 
scan in bias will need the high-energy features (such as 
the band curvature in the leads) to be well resolved as 
the contribution to the current from states with energies 
within the Fermi levels in the left and the right leads 
becomes important. Therefore, at AV ^ Tk, the best 
approach is to use A = 1 in the tDMRG calculations and 
subsequently perform a finite-size scaling analysis of the 
average currents. 

Our results are illustrated in Fig. [71 Fig.[7Ka) shows the 
current versus time for different bias values and A = 1, 
N ~ 72 and A = 2, = 32. As a general feature, the 
average current (J)t increases with the bias AV in all 
cases, also seen in Fig. [TJ^b). Moreover, it is evident 
that for the values of AV depicted in Fig. [ZKa), runs 
with either A = 1 or A = 2 give a qualitatively similar 
behavior. 

A more quantitative analysis of the {J)r{AV) curves 
obtained from the tDMRG data is presented in Figs. [ZKb) 
and[7{c). Deviations from perfect conductance ((J)r = 
GoAV) are clear at large bias (Fig. [7{b)). At small 
biases, the results are better visualized by numerically 
calculating the differential conductance d{J)r/d{AV), 
shown in Fig. [TKc). For small AV , this quantity is equiv- 
alent to the linear conductance. 

For A = 1 and N ^ 72 sites, we sec that the corre- 
sponding d{J)r/d{AV) curve is substantially below Gq 
for AV < Tk in sharp contrast with the A = 2 results, 
emphasizing the importance of using Wilson chains at 
small biases. At large biases, we have performed a finite- 
size scaling analysis with N < 80, and we conclude that 
the data displayed in Fig. [TJc) are converged down to 
AV > 0.15 with an uncertainty of S [d{J)r /d{AV)] - 
0.01. 

In an intermediate bias regime of Tk S ^ 0.15, 
the A = 1 results are still plagued by finite-size effects. 




0.1 0.2 0.3 0.4 

AV 
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FIG. 7: (color online) (a) Current J(t) versus time for 
AV = 0.1, 0.2, 0.3 and A = 1 (iV = 72) and A = 2 (A^ = 32). 
(b) Average current (J)t and conductance d{J).r /d{AV) as a 
function of bias. The dashed line is {J)t = GoAV. (c) Differ- 
ential conductance for A = 1 (A^ = 72) and A = 2 {N ^ 32). 
The dotted line represents an interpolation between the low 
bias regime (better described by A = 2) and the high bias one 
(better described by A = 1). 



while the A = 2 results overestimate the actual steady- 
state currents. In this regime, a discretization scheme 
that interpolates between linear and a logarithmic dis- 
cretization will likely be the best choice. Qualitatively, 
one can estimate the expected result by a linear interpo- 
lation between the A = 1 and A > 1 data, as illustrated 
by the dotted line in Fig. [TKc) 

Finally, note that the finite-bias regime of a single- 
impurity Anderson model has been studied with tDMRG 
by Kirino et al. fRef. [25l) for U /V < 2.8 and system sizes 
of up to A^ = 64 (all at A = 1). Here, we slightly exceed 
that regime by working at U/T ~ 3.125, and wc also con- 
sider larger system sizes of A^ = 64, 72, 80. Our key point, 
though, is that tDMRG runs at AV^ < Tk notoriously 
improve and correctly capture Kondo correlations when 
performed with Wilson chains . 



VII. SUMMARY 

In this paper, we applied the tDAIRG method to the 
study of transport through a quantum dot coupled to 
nonintcracting leads with a logarithmic discretization. 
This yields a considerable improvement over tDMRG 
studies with real-space tight-binding leads, as it extends 
the parameter space in which known exact results can 
be reproduced. One of the main advantages of the ap- 
proach is that smaller chains are sufficient to obtain the 
expected result of a perfect conductance for the single- 
impurity problem at particle-hole symmetry. 
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In spite of the challenges imposed by the longer time 
scales needed for the description of the Kondo regime, 
the study of transport properties in nanostructures with 
time-dependent DMRG brings several advantages over 
other methods: it is straightforward to adapt codes to 
more complicated geometries and time-dependent Hamil- 
tonians, including correlation effects in the leads, as well 
as systems far from equilibrium, such as transport be- 
yond the linear response regime (finite-bias) For the 
latter example, we presented a case study at an interme- 
diate U/V to argue that for bias values (AF < T^), a 
logarithmic discretization should preferentially be used, 
while at large bias, a regime in which high-energy fea- 
tures of the leads dominate, a tDMRG study with A = 1 
and a finite-size scaling analysis yields better results. 

Furthermore, we believe the general concept of utiliz- 
ing a logarithmic discretization in tDMRG can have a 
broader range of applicability in the description of Kondo 
systems. In particular, it can play a key role in the calcu- 
lation of nonequilibrium dynamic correlation functions, 
as recently highlighted in NRG-based approachesi^i 
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APPENDIX A: NONINTERACTING CASE 

In this appendix, we present exact diagonalization re- 
sults for the time-dependence of the current Eq. ([4]) for a 
chain of = 32 sites and t' = 0.4. Figure [5] shows J(r) 
vs. time for A = 1, 1.4, 2. Some main features induced by 
the discretization discussed in Sec. lIIII are already present 
in the noninteracting case: (i) as A is increased, the sign 
reversal of the current occurs at much later times and 
(ii) while the average current in the [/ = case remains 
at G = Go, the dispersion around this mean value is 
significantly enhanced by a A > 1. 

More specifically, by taking averages over suitable time 
intervals, we find G/Gq = (1.00 ± 0.04) and G/Go ^ 
(l.OOib 0.13), for A = 1.4 and 2, rcspectivlcy. The devia- 
tions are computed from Eq. We thus conclude that 



3 - u=0, t'=().4, N=32, ED, AV=1()'" 
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FIG. 8: (color online) Current J(r) as a function of time 
for the noninteracting case {U = 0, Vg = 0, t' = 0.4, A'^ = 
32) for different values of A = 1, 1.4, 2, calculated with exact 
diagonalization. 



the oscillations seen in Fig. [2] in the interacting case are 
due to the discretization. 

Moreover, from the noninteracting case we learn that 
using the logarithmic discretization, one manages to re- 
produce the average quite well. Once one has achieved 
that in the interacting case for a pair of ( A^, A) , wc expect 
that the oscillations will be reduced by increasing N and 
decreasing A at the same time in a controlled way such 
that the average current J(t) remains constant. 

An alternative way of suppressing the oscillations in- 
duced by the discretization is to exploit the so-called z- 
trickiiii^i^ Indeed, this works quite well: the thick solid 
line in Fig. [S]is obtained by using the z-trick for A = 2, 
which clearly improves the data quality over the simple 
A = 2 curve. However, it turns out to be necessary to 
average over many values of z: results shown in Fig. [5] 
were obtained by averaging over forty J(r)-curvcs with 
z- values ranging from to 1. In practice, this makes 
this procedure numerically expensive for tDMRG calcu- 
lations. 



APPENDIX B: COMPUTATIONAL ASPECTS 

As a key result of this work, we have argued that using 
the logarithmic discretization much smaller chains than 
in the A = 1 case can be used to obtain equally good, 
if not better, results for the conductance. This suggests 
a gain in the computational costs needed to obtain the 
numerical results, by reducing the required system sizes 
roughly by a factor 4. 

For a more stringent estimate of the computational ef- 
ficiency, we consider the entanglement growth during the 
time-evolution. We measure this quantity by computing 
the block entropy Si , associated with the reduced density 
matrix pi of a DMRG block of length I: 

Si^^{pMpi). (Bl) 
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FIG. 9: (color online) Block entropy Si for a block of length 
I = 16 vs. time for U jV = 5 (A = 2) and U jV = 8 (A = 
1, 2, 3). iV = 32 and AV = 10"^, the time steps (5r, and the 
truncated weight 5p are given in the legend. 

The reduced density matrix pi is obtained at each step 
by dividing the DMRG chain (the so-called super-block) 
into "system" (size /) and "environment" parts, and trac- 
ing out the environment's degrees of freedom (see, e..g.; 
Ref. [3^ for a discussion of the DMRG method) . An in- 
crease in 5"/ renders a simulation inefficient as time or 
system size grows, since more DMRG states need to be 
kept in order to keep the truncation error below a given 
threshold. 

In Fig. [HI we plot the block entropy for a block of length 
Z = 16 as a function of time, for two different values of 
C//r = 5 (A = 1) and 8 (A = 1,2,3). Typically, the 
entropy rapidly increases at short times, but at times r > 
3, it exhibits a linear increase in time for A > 1, S"/ oc r. 
This is the expected behavior for a global quench^^i^iS 
yet it is a non-obvious one here as the excitations in the 
leads do not travel at a constant velocity (see Fig.[3I^b)). 
The oscillations in Si(t) seen in the A = 1 case (dotted 
line) are due to the sign reversal of the current. The 
key point is that this linear increase is slow, i.e., the 
prefactor is small. During the time interval r € [5, 100], 
Si only grows by a few per cent. This is ultimately the 



reason why we can push our tDMRG runs to times long 
enough to reach the steady state for U jV < 6 at moderate 
numerical costs, especially since the entanglement growth 
is the weaker the larger IJ jV is. 

We thus observe that the entanglement growth, i.e., 
the increase in the entropy S*;, depends on both A and 
U jV . Yet, it is fortunate that in the case where longer 
times are needed in order to capture the steady-state 
current (large IJ /F) the increase in 5*; is weaker. 

It is illustrative to give an example on what the en- 
tanglement growth implies in practice for the numerical 
effort when working at a fixed truncation error Ep. For 
C//F = 3.125, we find it sufficient to keep to ~ 280 states 
at A = 1 in order to ensure a maximum truncated weight 
of &p ^ 10^^ on a chain of = 32 sites, compared to 
TO ~ 1000 at A = V2 and TO ~ 1600 at A = 3 (both num- 
bers refer to times t < 30 with = 0.05). At a larger 
L//F, say 12.5, this relaxes to to '--^ 200 and to ^ 400, for 
A = 1 and A = -\/2, respectively. 

We finally comment on the generic tDMRG errors, 
the accumulated truncation error, and the Trotter error. 
These are not independent: the smaller the time step, 
the faster truncation errors will accumulate^S We justify 
our choice of parameters by considering the numerically 
worst case, i.e., a small J7/F ~ 3. 

For most calculations, keeping a maximum of to = 660 
states up to times of order r ^ 50 is sufficient to keep 
the truncation error below 10~^ (^ 3 x 10~^ in a few 
sweeps). More importantly, we have checked that, keep- 
ing up to TO = 1600 states, the current J(t) is practically 
converged: for the case of Fig. [^fc) and A = 2, the max- 
imum relative change in J(t) is ^ 1%, comparing runs 
with 8p = 10^^ and 8p = 10~^. In addition, we have 
calculated the forth-back error— to validate that this is 
a sufficiently small discarded weight for our purposes, 
in which the oscillations cause the dominant fluctuation 
around the current's average (see App. 

We have further checked our tDMRG with the chosen 
parameters against exact diagonalization for the nonin- 
teracting case to make sure that the so-called ruir-away 
time^ is not the limiting factor in our case. 
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